function uh = hybrid_Limiter(uh)

global Nx dimPk bcL bcR hx phiG NumGLP
uhb = [[0,0,0];uh;[0,0,0]];
uhmod = zeros(Nx,dimPk);
uhmod(:,1) = uh(:,1);
eta=zeros(Nx);
beta=1;

uhG = zeros(Nx,NumGLP);

for i = 1:Nx
    for d = 1:dimPk
        uhG(i,:) = uhG(i,:) + uh(i,d)*phiG(:,d)';
    end
end

% set_bc
if bcL == 1
    uhb(1,:) = uh(end,:);
end

if bcR == 1
    uhb(end,:) = uh(1,:);
end

for i = 1:Nx
    deltaUR = uh(i,2) + (2/3)*uh(i,3);
    deltaUL = uh(i,2) - (2/3)*uh(i,3);
    deltaURM = uhb(i + 2,1) - uhb(i + 1,1);
    deltaULM = uhb(i + 1,1) - uhb(i,1);
    
    deltaURM1 = minmod(deltaUR,deltaURM,deltaULM);
    deltaULM1 = minmod(deltaUL,deltaURM,deltaULM);
    
    uhmod(i,2) = (deltaURM1 + deltaULM1)/2;
    uhmod(i,3) = 3*(deltaURM1 - deltaULM1)/4;
 end

%     for i=2:Nx
%              eta(i)=eta(i)+abs(uhG(i,3)-uhG(i-1,3));  
%     end
%         for j=1:dimPk
%              eta(i)=eta(i)+abs(uhG(1,3)-uhG(Nx,3));
%         end 


    for i=2:Nx
        for j=1:dimPk
             eta(i)=eta(i)+abs(uh(i,j)-uh(i-1,j));  
        end
    end
        for j=1:dimPk
             eta(i)=eta(i)+abs(uh(1,j)-uh(Nx,j));
        end 


  for i=1:Nx
        eta(i)=eta(i)/(2*hx*beta);  
  end
       Max=max(eta(:));

  for i=1:Nx
      if(eta(i)<1)
      else
      uh(i,:)=(Max-eta(i))/(Max-1)*uh(i,:)+(1-(Max-eta(i))/(Max-1))*uhmod(i,:);
      end
   
   end

end

function a1 = minmod(a,b,c)

global hx M

if abs(a) < M*hx^2
    a1 = a;
else
    if sign(a) == sign(b) && sign(a) == sign(c)
        a1 = sign(a)*min(abs([a,b,c]));
    else
        a1 = 0;
    end
end
end




